Tailoring laser pulses with spectral and fluence 
constraints using optimal control theory 



J. Werschnik and E.K.U. Gross 

Institut fiir Theoretische Physik, Freie Universitat Berlin, Arnimallee 14, 14195 
Berlin, Germany 



in 
o 
o 

(N 

E-mail: j an . werschnikOphysik . f u-berlin . de 

<' 

(N 



Abstract. Within the framework of optimal control theory we develop a simple 
iterative scheme to determine optimal laser pulses with spectral and fluence constraints. 
The algorithm is applied to a one-dimensional asymmetric double well where the 
, control target is to transfer a particle from the ground state, located in the left well, 

to the first excited state, located in the right well. Extremely high occupations of the 
first excited state are obtained for a variety of spectral and/or energetic constraints. 
I Even for the extreme case where no resonance frequency is allowed in the pulse the 

■ algorithm achieves an occupation of almost 100%. 

o 

:^ 
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, , 1. Introduction 

■ In the last fifteen years, tlie control of quantum meclianical systems via liglit pulses lias 
experienced/seen considerable progress, both on the experimental and on the theoretical 
side. Quantum control experiments have been pushed forward by the improvement 
of laser pulse shaping jH 13 E] and the implementation of closed-loop learning (CLL) 
techniques Experiments using CLL delivered highly encouraging results, ranging 
from the control of chemical reactions El 13 13 UHl UH HI] to the control of high- 
harmonic generation ITi] . 

On the theoretical side, the most important contributions have been the 
introduction of optimal-control theory [13 CEl HZI and the continuous development 
of rapidly converging iteration schemes [IHl UHl 1201 to calculate optimal laser pulses. 
Recently, some of these schemes have been generalized to include dissipation ^21j, to 
account for multiple objectives [22] and to deal with time-dependent control targets 

[23,21123]. 

Most fruitful are investigations where theory and experiment come together: 
The theoretical analysis of laser pulses can be useful, and sometimes is essential. 
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in deciphering the pulse shapes obtained from experiment To speed up the 

convergence process of the experimental learning cycle, calculated pulses can be used 
either to provide an initial guess or to reduce the gigantic search space by determining 
the most important shape parameters j2Zl • Besides these direct applications, theory and 
computer simulations make it possible to explore the feasibility of future experiments 
and help to determine the requirements on the laser system and the pulse shaping 
device. Computer simulations can also help in understanding new ultrafast transition 
processes in laser-assisted chemistry |28[ EH! and in developing new implementations for 
the quantum computer jSHl ISI] • 

For all these applications it is extremely important that the computational schemes 
are able to include experimental constraints, such as limitations on the spectral 
bandwidth and on the laser fluence, i.e. the time-integrated intensity. As discussed 
in [Appendix A.l} a pulse from an unconstrained optimization will perform much worse, 
if the constraint is applied "brute force" after the optimization than a pulse coming 
from a scheme where the same constraints are built in. 

So far, only few attempts have been made to take restrictions of this kind into 
account. In reference a scheme to calculate the pulse for a given fluence is shown. 
However, it does not make use of the immediate feedback introduced in [TH] and suffers 
from a rather unstable convergence. A constraint on the spectrum is considered in 
|32] for a steepest descent method which, in the quantum control context, is found to 
suffer from poor convergence and a strong dependence on the initial pulse jHSl- An 
elegant way to restrict the spectrum has been presented by the authors of reference [Blj. 
This scheme preserves the rapid and monotonic convergence behavior of the underlying 
scheme [IB] by projecting out undesired parts of the time-dependent wave-function, 
which are responsible for the unwanted spectral components. However, this method is 
not sufficiently general and does not easily allow for an additional fluence constraint (as 
it keeps the unphysical penalty factor). 

The scheme presented in the following allows one to incorporate fluence and/or 
spectral constraints in the optimization and it shows very good convergence, when 
applied to a ID model, although a proof of monotonic convergence similar to reference 
[To] does not go through here. Further more, the scheme is very simple to implement. 

An introduction to optimal control theory is given in section |21 We then explain 
our schemes in section |3] In section HI we present a test system and discuss the 
numerical details. The results from applying our algorithms to this system are analyzed 
in section |31 

2. Optimal Control Theory 

In this section we sketch the basics of optimal control theory applied to quantum 
mechanics. We consider an electron in an external potential ^(r) under the influence 
of a laser field propagating in z-direction. Given an initial state \I^(r, 0) = 0(r) the time 
evolution of the electron is described by the time-dependent Schrodinger equation with 
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the laser field modeled in dipole approximation (length gauge) 

z^vl>(r,t) = #vl/(r,t), (1) 

H =Ho-fie{t), (2) 

^0 =f + V, (3) 

(atomic units are used throughout: h = m = e = 1). Here, fj, = {jix,jly) is the dipole 



operator, and e{t) = (exit), ey{t)) is the time-dependent electric field. The kinetic energy 
operator is T = — 

Our goal is to control the time evolution of the electron by the external field in such a 
way that the expectation value of the target operator O is maximized with respect to 
the wave function at the end of the pulse \l/(r, T). Mathematically, this goal corresponds 
to maximizing the functional fSl QHl : 

jm = {^iT)\d\<f{T)). (4) 

Usually O is assumed to be positive-semidefinite which guarantees monotonic 
convergence of the schemes discussed in references |18| IT?I| EDI 120] • A few examples 
will be discussed at the end of this section. 

The functional Ji[^] will be maximized subject to a number of physical constraints. The 
idea is to cast also these constraints in a suitable functional form and then calculate 
the total variation. Subsequently, we set the total variation to zero and find a set of 
coupled partial differential equations [131111 • The solution of these equations will yield 
the desired laser field e(t). 

In more detail: Optimizing Ji may possibly lead to fields with very high, or even 
infinite energy. In order to avoid these strong fields, we include an additional term 
in the functional which penalizes the fiuence of the field. This can be done for each 
polarization direction separately: 

rT 

M(^]= - / dtaje/{t) j = x,y, (5) 

where aj is a penalty factor that has to be chosen. It balances the optimization between 
increasing the yield and restricting the energy to achieve the maximal value for the com- 
bined functional Ji + J2. Note, the penalty factor aj can be made time-dependent to 
restrict the laser pulse to a certain shape 



The constraint on the laser fiuence can be expressed also in another way: 



E 



dt e/(t) - E^ 



0, 



(6) 



Here, aj is a (time-independent) Lagrange multiplier. Instead of specifying aj we have 
to prescribe specific values, Eq., for the components Eq^ and Eq^ of the laser fiuence. 
Hence, this approach requires two Lagrange multipliers and ay. 
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The constraint that the electronic wave function has to fulfill the time-dependent 
Schrodinger equation is expressed by 

Me, ^,x]= - 253 dt (xW I [idt - if) I ^{t)^ . (7) 

with a Lagrange multiplier x(r,t). \E'(r, t) is the wave function driven by the laser field 
e{t). 

The Lagrange functional has the form 

J[X, ^, e] = + Me] + Mx, e]. (8) 

Setting the variations of the functional with respect to X; and e independently to 
zero yields 

afyit) = - Q{x{t)\M'^{t)), j=x,y (9) 
= («9t-^)^(r,t), ^(r,O) = 0(r), (10) 

(idt-H^ X{r,t) = 

z(x(r,t)-6vi>(r,t))5(t-T). (11) 

Equation © determines the field from the wave function \E'(r, t) and the Lagrange 
multiplier x(r, ^)- 

Equation (fTUI) is a time-dependent Schrodinger equation for \I'(r, t) starting from a given 
initial state 0(r) and driven by the field e(t). If we require the Lagrange multiplier x(r, t) 
to be continuous, we can solve the following two equations instead of PHI : 



Wt-Hjx{r,t) = 0, (12) 
x(r,T) =OM/(r,T), (13) 

To show this we integrate over (|T0|) 

lim [ dt \(idt - h) x(r,t) 

= lim / dti(x{r,t)-d{t)m{Y,t)]6{t-T). (14) 

The left-hand side of (I14j) vanishes because the integrand is a continuous function. It 
follows that also the right-hand side must vanish, which implies fll3|) . From equations 
ffT8|l and (HU then follows equation 

Hence, the Lagrange multiplier x{^ii) satisfies a time-dependent Schrodinger equation 
with an initial condition at t = T. The set of equations that we need to solve is now 
complete: ©, (fTnjl . (fT^ and (fT!?|) . If we use J2 instead of J2 we also have to perform a 
variation with respect to Uj which simply yields the restriction: 

I dte%t) = Eo^. (15) 
^0 
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To find an optimal field e(t) from these equations we use an iterative algorithm which 
is discussed in the next section. 

We conclude this section with a discussion of the target operator O. Basically, 
there exist two classes of target operators. Namely, operators that are non-local, 
e.g. projection operators and operators that are local (multiplicative), like the density 
operator. If we want to maximize the occupation of a given target state at the end 
of the laser pulse, we choose jTHl EE] , 

0=|$/)($/|. (16) 

This scheme can be extended to achieve multiple goals, i.e. to have different states 
populated at the end of the pulse. In that case one uses, 

d = J2(^k\¥>k){¥>k\. (17) 

k 

The factors (3k allow for the possibility to "fine-tune" the target occupations among 
each other (multi-objective optimization), i.e. to balance between the importance of the 
individual targets. For example, if we choose /?„ negative the optimization will avoid 
the occupation of the state |v?n)- 

Note, that this kind of multi-objective optimization is different from the target 
to reproduce a (coherent) superposition of field free eigenstates of the Hamiltonian Hq 
given by 

d =|$/)($/|, 

1$^)= X^Cfclv^fc). (18) 

k 

Using as target operator the projection operator (fTHI) leaves the freedom of a purely 
time-dependent phase factor for the wave function \E'(T). It is possible to fix the phase 
with the following functional: 

min ll^(T) - $jf = 2 (1 - 3ft(*(T)|$j)) (19) 

^ Ji = max3?(^(r)|$j), (20) 

where we have assumed normalization: (\I'(r)|\I'(T)) = = 1. 

The target operator may also be local jTHI- If we choose O = 6{r — fq) (the density 
operator), we maximize the probability density in fq at t = T: 

J^ = {^iT)\^\^iT))=n{ro,T). (21) 

For this control target, the optimization process will try to concentrate the density in the 
point tq at the end of the pulse |HE]- Numerically, the (5-function can be approximated 
by a sharp Gaussian function. 



3. Algorithm 



In this section we present iterative schemes for the optimization of laser fields under 
additional constraints on the fluence and/or on the spectral distribution. 
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3.1. Fluence constraint 



We first describe the algorithm which yields an optimized laser pulse producing an 
assigned value of Eq- (for each polarization direction j = x,y, cf. equation The 
set of coupled equations to be solved is given by equations Q, (fTUj) . (fT^ and (fTSj) . The 
scheme below shows the order in which these equations are solved in the kth iterative 
step. 



k-th step: ^('=)(0) 



[^(fc)(T) 



^(fc)(0)] 



(22) 



with the laser fields e^^\t), e''''\t) given by 

1 



j 



it) 



a 



{X^'\t)\fi,\¥'\t)), 



it) 



a 



a 



(fc) 
(fc+i) i 



it), 



J =x,y, 



(23) 



(24) 



where the Lagrange multiplier a^'^^^'' is defined by: 



a 



(fc+i) 



lodt 



a 



(kMk) 



it) 



(25) 



The initial conditions in every iteration step are 

^(r,O)=0(r), (26) 

x(r,T) = Ovl/(r,T). (27) 

The scheme starts with the propagation of \&^'^)(r, t) forward in time using the laser 
field e*^^-' (t) which has to be guessed. The result of the propagation is the wave-function 
\I^(°)(r, T) which is now used to calculate x^^K''^, T) by applying the target operator (j2Zj). 
We continue with propagating x'-"-* (r, t) backwards in time using the laser field 'e^^^t) 
fl23|) . To solve equation (j^Hj) we have to know both wave functions \l/'^°)(r, t) and x^'^\'^, t) 
at the same time t, which makes it necessary to either store the whole time-dependent 
wave function \l/(°)(r,t) or propagate it backwards with the previous laser field e'^°)(t). 
The version of the algorithm that avoids storage is indicated by the brackets in the 



scheme 



Besides that, it is necessary to provide an initial value for of* which we 



choose to be: 



a 



(0) 



Kdt ef\t) 



En 



The result of the backward propagation x^^K'^, t) is the laser field e^^\t) which we now 
re-scale to the right value yielding e^^\t). This completes the first step. The second 
{k = 1) or, in general, the kth iteration repeats the described procedure starting again 
from the initial state \l/*^'^)(r, t) = 0(r) and applying the rescaled field e^'^^t). 
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The scheme described above has some aspects in common with the techniques 
described in reference ^3] and reference ^H]: The basic idea of incorporating fluence 
constraints in the optimization algorithm was given in reference j.l5]. However, the 
authors do not make use of the immediate feedback (cf. equation ((221)), i-e. the backward 
propagation is accomplished by updating x(r, t) and e{t) in a self-consistent way, which 
was suggested in reference [18j. On the other hand, the technique presented by the 
authors of reference [TH| does not allow to build in fluence constraints, since aj is not a 
Lagrange multiplier in their case. Roughly speaking, the technique presented above is 
a combination of both approaches. 



3.2. Spectral constraint 

The algorithm with built in spectral restrictions is similar to the one presented above 
with two important differences: The factor aj is a penalty factor. It has to be specified 
from the start and remains unchanged during the optimization. Second, the update of 
the field e^'^'''^'' (t) in equation pH) has to be replaced by: 



ef^'\t) = T X T f;\t)\\ J = x,y, (28) 

where the symbol JF indicates a Fourier-transform. The spectral constraint is formulated 
in terms of a filter function fj{uj). Since ej{t) is real valued we have to make sure that 
fj{uj) = fj{—u). For example, the filter function could be chosen to be: 

fj{uj) = exp[-7(u; - uoY] + exp[-7(u; + ujo)^], (29) 

so that only the frequency components around the center, icuo, of the Gaussians are 
allowed in the pulse. If one uses instead: 

fj{u) = 1 - (exp[-7(a; - uq)^] + exp[-7(u; + ujof]) , (30) 

one would allow every spectral component in the laser field except the components 
around ±uJo. 



3. 3. Spectral and fluence constraint 



Finally, we note that both schemes can be combined. This combination makes it 
possible to incorporate even more realistic experimental constraints in computational 
pulse optimizations. This is achieved by the scheme (j22|l . the equation p3| and: 



fit) 



a 



a 



j,{^)T^\t) 



J =x,y, 



{k+i) 



(31) 
(32) 



where af~^^^ is evaluated with the filtered field e'j^\t) 



a 



(fe+i) 



ffdt 



(k)-{k) 

C(- €■ 
3 3 



En 



(33) 
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to yield the right value for Eq.. The total spectral power is related to the fluence 
Parseval's theorem: 

/ + 00 1 Z' + OO 

dtemT-t)[e,it)f = 7r / 
■oo J -oo 

In this combined form we first apply the filter function to the laser field (|HT|) then 
we rescale the field to yield the right value for Eq. (see equation (jHH))- 
We conclude the section with a few remarks: 

• For each polarization direction one can specify a separate filter or fluence. 

• The convergence proofs of references fHl ^] do not go through in our case. This 
is due to the changing value for aj'^ and, in the case of spectral constraints, due to 
the modified field ()28j) . However, as will be shown in section EJ we still find a very 
good convergence of the presented algorithms in the numerical examples. Even for 
the "brute-force" spectral filter we find a satisfying convergence behavior (unless 
too many essential features of the pulse are suppressed by the function fj{uj)). 

• Since we do not expect a monotonic convergence we have to add some additional 
intelligence to the algorithm, e.g. we store the field which produces the pulse with 
the highest yield and consider this field as the result of the optimization. 

4. Computational details and model system 

We choose a one- dimensional asymmetric double well to test our algorithms. The double 
well is similar to reference but has an additional cubic term: 

V{x) = ^x'-'^x' + Px\ (34) 

with ujq corresponding to the classical frequency at the bottom of the well and the 
parameter B adjusting the barrier height. The number of pairs of states below the 
barrier is approximately B/ujq. Here, we choose 5 = cuq = 1.0 and j3 = 1/256 which 
leads to two states below the barrier, as shown in figure H In order to analyze the laser 
pulses from the optimization runs we calculate the excitation energies (see table and 
dipole moments (see table ^ of the system by propagating in imaginary time. 

Table 1. Excitation energies in atomic units [a.u.] for the ID asymmetric double well, 
calculated by imaginary time propagation. 





|0) 


|1> 


|2) 


|3) 


|0) 


0. 








ll> 


0.1568 


0. 






|2) 


0.7022 


0.5454 


0. 




|3) 


1.0147 


0.8580 


0.3125 


0. 


|4) 


1.5294 


1.3726 


0.8273 


0.5147 



Spectral and fluence constraints using optimal control theory 



9 



0.5 



7 

S -0.5 



-1 

-6 -4 -2 2 4 6 
X [a.u.] 

Figure 1. The plot shows the model potential with the ground-state ( ), the 

first excited state ( ), the second excited state ( ) and the third excited state 

( — • — ). Each state is shifted according to its eigenvalue. 

Table 2. Dipole matrix elements for the ID asymmetric double well, calculated by 
imaginary time propagation. 





|0) 


|1> 


|2> 


|3> 


|4> 


|0> 


-2.5676 










ll) 


0.3921 


2.3242 








|2) 


0.6382 


-0.7037 


-0.5988 






|3) 


-0.3865 


-0.4630 


1.7051 


0.1958 




|4) 


-0.1414 


0.2118 


0.1593 


-1.7862 


-0.0939 




The time-dependent Schrodinger equation for the ID double well is solved on a 
grid, where the infinitesimal time-evolution operator is approximated by the 2nd-order 
split-operator (SPO) technique [HH] : 

/ i-t+At 

fti exp(-^f At)exp(-z1/(t) At) x exp(-^ f At) + O(At^). 

Following the scheme described in section 01 one needs three propagations per iteration 
(if we want to avoid storing the wave function). Within the 2nd order split-operator 
scheme each time step requires 4 Fast Fourier Transforms (FFT) [SH] for the backward 
propagations, because we have to know the wave-function and the Lagrange multiplier 
in real space at every time-step to be able to evaluate the field from equation ©. For 
the forward propagation we only need 2 FFTs. This sums up to 10 FFTs per time step 
and iteration. 

The parameters used in the runs are summarized in table El The initial guess for the 
laser field was e^°''(t) = —0.2 in all calculations. This choice is arbitrary but has the 
advantage of producing a significant occupation in the target state at the end of the 
pulse, necessary to get the iteration working. Although the simple choice e*^°''(t) = 0.0 
will work as well in most cases, it represents a minimum of the functional since initial 
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Table 3. In ordinary runs the listed numerical parameters (given in atomic units) were 
employed. For the scans we have used a coarser grid in space and time, as indicated 
in the second column. 



parameter 


single run 


scan 




T 


400.0 


400.0 


pulse length 


•^max 


30.0 


20.0 


grid size 


dx 


0.1172 


0.1563 


grid spacing 


dt 


0.001 


0.005 


time step 


e(0) 


-0.2 


-0.2 


initial guess 



and target state are orthonormal. Therefore the algorithm could get stuck in principle. 
The obtained solutions, which are presented in the following chapter, are all far away 
from the initial guess. This suggests that the solutions do not depend on the initial 
guess for the laser field. 



5. Results 



In this section we apply the algorithms described above to our ID model for electron 
transfer. We start in the ground state |0) {t = 0) where the electron is localized in the 
left well and demand that at the end of the laser pulse {t = T) it will be transfered to 
1st excited state |1), which is mainly located in the right well (see figure^. The target 
operator in this case is a projection operator onto the first excited state: O = |1)(1|. 
Therefore, the success is measured by | (^I/(T) 1 1) p which we simply refer to as the "yield" . 
The pulse length is chosen to be T = 400 (^ 9.7 fs). 



5.1. Fluence constraints 

5.1.1. Fixed fluence. In the following we first apply our algorithm to find an optimal 
field with the fluence Eq = 0.080. This is the value obtained by an estimate using 
the two- level system (see Appendix A. 2D . After 894 iterations we obtain a yield of 



99.91% which is higher than the yield found by the two-level estimate. The optimal 
laser field and its spectrum (Fourier transform) are shown in fig ures |2(I)1 and |2(bj] 
The spectrum is dominated by three narrow peaks which correspond to the excitation 
energies Uqi = 0.156, uji2 = 0.545 and ujq2 = 0.702. This suggests that the optimized 
transition process is a mixture of the direct process, i.e. the excitation from |0) — > |1) 
and an indirect process which uses the second excited state as intermediate state: 
|0) 1 2) |1). Other indirect processes, like |0) — >• |3) |1), play only a minor 
role in this case. This interpretation is supported by looking at the evolution of the 
occupation numbers in time (figure |2(c)] ) . First, the laser pulse populates the second 

excited state ( ) and then after half of the pulse duration depopulates it again. 

Looking once more at the spectrum (in figure |2(bj| ) we observe a group of peaks around 
Uqi (t^ £ [0,0.1] and u G [0.2,0.3]) which do not correspond to any excitation energy 
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of the field-free Hamiltonian. However, these frequencies play an important role in 
the transition process. If we filter out these frequency components, rescale the fluence 
to Eq = 0.080, and then propagate this modified laser pulse, we find, at the end of 
the pulse, the following occupations: ground-state 16%, first excited state 38%, second 
excited state: 44%, and in all higher levels 6%. In particular, the direct transition 
and the back transfer from the intermediate level |2) to the target state in the indirect 
process are less efficient without these extra frequencies. Further analysis of this kind 
shows that the low-frequency components and especially the zero frequency component 
(bias) are crucial since they introduce a (slight) shift of the resonance frequencies, visible 
as a broadening of the uooi peak in figure |2(b)[ If, on the other hand, these components 
are missing the remaining frequencies become slightly off-resonant, resulting in the low 
efficiency of 44%. 

If we filter out everything except the extra peaks we find a target state occupation 
of 1%. Understanding these extra peaks as a third type of transfer process (see 
section I5.3.3j) suggests that, in this case, a mixing of transition processes seems to 
be superior in terms of the maximum target yield than a pulse consisting of a single 
process only, e.g. the direct process. 

The final yield 99.91% is only 0.61% better than the yield coming from the 
simple monocromatic pulse estimate of the two-level system. This gain has a high 
price, the optimized pulse is hardly realizable in any experiment. Although the gain 
improves with shorter pulse lengths (see Appendix A.2 ), this example demonstrates the 



typical dilemma between theory and experiment: Calculated pulses often have a far too 
complicated spectrum to be produced in practice. In section \E7]\ and we demonstrate 
how this dilemma can be resolved. 

To conclude the analysis we look at the convergence behavior of the applied scheme 
(see figure gdj. We find a fast convergence within the first 20 iterations. After these 
20 iterations the improvement of the yield slows down, like it is also found in the rapid 
monotonic schemes presented in references [THl 120] • 

5.1.2. Energy versus yield We apply our method to scan through a range of values for 
Eq from 0.010 . . . 1.000. The scan, displayed in figure El shows that there seems to be 
a critical value Eq which is necessary to get very high occupations (|(\Ef(T)|l)p > 0.99 
of the target state. For values Eq > Eq the algorithm always finds a laser field that 
produces yields above 99%. 

For long pulse durations (as it is the case here) we can give a rough estimate of this 
critical value Eq with the help of the two-level system: 

^"'=4 = air (35' 

(36) 

with A = Ti/ (fiQiT) (see [Appendix A~ ). 

If we take a closer look at some of the optimized fields (see figure 4(a)| ) for the 
values Eq = 0.010,0.050,0.100,0.200, we see that the spectra (see figure |4(bj| ) of these 
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Figure 2. We apply the algorithm for Eq = 0.080. The optimized field is shown 
in (a), with its spectrum in (b). In (c) we plot the time evolution of the occupation 

numbers (n = ( ), n = 1 ( ),n = 2 ( ) and n = 3 ( )). 

The convergence behavior is shown in (d). Also shown in (d) are the values of the 

Lagrange multiplier a ( ) during the iteration which we scaled by a factor of 0.1. 

The O indicates the iteration with the highest yield. 



pulses get more complicated as the assigned fluence increases. In the lower two panels 
{Eo = 0.010, 0.050) we flnd peaks at the exact resonance frequencies. The optimized 
fields result in occupations of 35.24% and 97.63%. While the pulses shown in the two 
upper panels {Eq = 0.100,0.200) produce yields of 99.81% and 99.95%. The peaks 
corresponding to the direct |0) — > |1) and indirect process |0) — > |2) |1) are "Stark" 
shifted. For the stronger pulses, we also find an increasing low-frequency part. The 
spectrum in the top panel {Eq = 0.200) is difficult to analyze. However, one can see 
that more and more processes are taking part in the transition, i.e. peaks occur near 
the other resonance frequencies. 



5.2. Spectral constraints 

In the following we present the results of the algorithm with spectral constraints and 
penalty factor for two examples of the filter function. These examples are motivated by 
the findings of the previous chapter, namely that the transfer of the particle occurred 
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Figure 3. The graph shows the yield for different values of Eq ranging from 0.010 
to 1.000. Each point corresponds to a single optimization. For these rmis we have 
used a smaller grid {xmax = ^Xmin = 20.0, dx = 0.15625) and a larger time-step 

{dt — 0.005). The vertical line ( ) corresponds to Eq — 0.080. Beyond this line 

only yields higher than 0.99 are found. 
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Figure 4. In (a), we plot the optimized laser fields for different values of Eq = 
0.010,0.050,0.100,0.200. Graph (b) shows the corresponding spectra. The resonance 
frequencies for the transitions |0) |1), |1) |2) and |0) |2) are indicated by 
vertical lines. In the two upper plots of graph (b) the peaks are Stark-shifted. Note, 
that these spectra also contain a large low frequency part. 



via a mixture of a direct transition and indirect transitions. We want to find a laser 
pulse that produces a high yield and only contains spectral components centered around 
the resonance frequency uqi. We know that such a pulse exists, since it appeared in 
the second iteration when looking for a pulse with Eq = 0.080 (see section IB.l.lj) . In 
the second example we optimize a laser that is not allowed to contain the excitation 
frequency uqi of the direct process. 



5.2.1. Direct transition Using spectral constraints in the optimization scheme allows us 
to explicitly select the direct transition, i.e. we search for a pulse whose main frequency 
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component is the excitation energy cuqi- This is done by applying a Gaussian shaped 
frequency fiher /(tu), according to equation (j^ . centered around = tuoi and with 
7 = 500. 

After 50 iterations the algorithm finds a laser pulse which results in a yield of 
99.97% . We set the penalty factor a = 0.05 and obtain a value of Eq = 0.090 which 
is slightly higher than the estimate from the two-level model but also more effective. 
The slight envelope on the field, shown in figure |5(a)[ stems from the finite width of 
the Gaussian (see figure |5(bj| ) . Frequency components near uqi are still allowed in the 
pulse and result in a beat pattern. The time dependent occupation numbers confirm 
that the higher states are not occupied during the transition (see figure |5(cj| ). The 
convergence, shown in figure |5(cj| is rather smooth. Note, that if we desire a sinusoidal 
field with a constant envelope we have to reduce the width of the Gaussian to allow 
only one single component in the spectrum (a Kronecker delta). Using such a filter 
we obtain a yield of 99.79% and Eq = 0.085. The field oscillates with the amphtude 
A = 0.0207 which is slightly higher than the amplitude derived from the two-level system 
Appendix A.'2). 



see 



5.2.2. Forbidden direct transition By choosing the complement of the filter function 
from the previous example, i.e., by allowing every frequency component except c^oi, we 
can optimize a field which also produces a very high yield. The filter function is given 
by equation (jHUjl with ujq = ujqi and 7 = 500. 

The optimization procedure (with a penalty factor a = 2.5) results in a target 
state occupation of 99.60% after 269 iterations. The optimized laser field is presented in 
figure |6(a)| it integrates to a fluence of Eq = 0.130. Its spectrum, shown in figure |6(b)| 
consists of two major components: Ua = 0.581 and = 0.676 which correspond to 
the Stark-shifted excitation energies uju and ujq2, i.e. the optimization takes care of the 
frequency shifts introduced by the large bias (zero-frequency component) of the field. 
That the transition occurs via the indirect process is confirmed by looking at the time- 
dependent occupation numbers, shown in figure |6(cjj First, the field starts populating 
the second excited state and then transfers the population to the target state. Other 
indirect processes, e.g. |0) — > |3) |1) or |0) — ^ |2) |3) |1), play only a minor 
role: The occupation of the third excited state stays below 2.5% and the frequency 
components correponding to these processes are very small. 



5.3. Combination of Spectral and fluence constraints 

The next examples demonstrate that even more restrictions are possible and we can still 
obtain very good yields. We combine the spectral restriction with the fluence constraint 
and continue the above examples by selecting among the indirect processes. Only two 
frequencies are allowed in the laser pulse and in addition we fix the fluence. In the last 
example we show that it is not even necessary to have resonance frequencies inside the 
laser pulse to reach very high occupations of the target state. 
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Figure 5. We apply the optimization algorithm but allowing only a Gaussian 
frequency distribution around luqi — 0.1568. The resulting field is displayed in graph 

(a). The filter function /(w) ( ) which is scaled by 0.01 is shown together with the 

spectrum in (b). The time dependent occupation numbers (c) confirm that only the 

ground state ( ) and the first excited state ( ) take part in the process. The 

second excited state population ( ) is hardly visible. The convergence is shown in 

graph (d). The O indicates the iteration with the highest yield. 



5.3.1. Selective transfer via intermediate state |2) In the former examples we found 
that the indirect process |0) — > |2) — > |1) plays a major role in the excitation process. 
Since it appeared always together with other processes, e.g. in section 15.1.11 together 
with the direct process or in section IB . 2 . 21 together with other indirect processes, we try 
to find a laser field with only the two excitation energies ujq2 and uJi2 (7 = 500) and in 
addition require Eq = 0.160. For these high requirements we have to pay a price which 
is the irregular behaviour of the yield during the iteration, shown in figure |7(d)| After 
540 iterations we find a yield of 99.90%. The restriction of the laser frequencies results 
exactly in the desired transition process, which is confirmed by the time-dependent 
occupation numbers, shown in figure [7(c)] 



5.3.2. Selective transfer via intermediate state |3) The process |0) — |3) |1) using 
the third excited state as intermediate state played only a minor role in the examples 
considered so far. Here, we try to optimize the laser pulse so that the transition is only 
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Figure 6. We prohibit the direct transition from |0) |1) by using the complement 
of the previous fiher function /(lu) = 1 — /{lo). The optimized field is shown in graph 

(a). The filter function f{uj) ( ), scaled by 0.01, is plotted together with the 

spectrum in (b). The time-dependent occupation numbers (c) confirm that now the 

second excited state ( ) plays a major role in the transition (ground state ( ), 

first-excited state ( ) and third excited state ( — • — )). The convergence is shown 

in graph (d). The O indicates the iteration with largest occupation of the target-state. 

performed via this process. In addition we require Eq = 0.320. Again, we use a double 
Gaussian filter, one Gaussian centered at cuis, the other one at cuqs choose the width 
parameter 7 = 500. 

The results are shown in figure |H1 Like in the previous example, the high 
requirements on the laser field result in a rather erratic convergence (see figure |8(d)] ) . 
The field, shown in figure |8(a)[ produces a target state occupation of 99.89% after 162 
iterations. The time-dependent occupation numbers (see figure |8(c)| ) show that the 
transition exactly happens in the desired way. 

5.3.3. Low frequency pulse Using a band filter, i.e. f{uj) = 6{uj — Ua) * 0{ujb — 
u) + 9{—uJa — uj) * 6{ijj + tUfe), with uja = 0.000 and = 0.120 we can find a laser 
pulse resulting in high yields without allowing any resonance frequency in the laser 
spectrum. The allowed frequencies {u G [0.000,0.120]) are smaller than the lowest 
excitation frequency Uqi. The additional constraint on the fluence in the optimization 
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Figure 7. Here we use a double gaussian window (7 = 500), allowing only frequencies 
around a;o2 and LU21 and in addition demand Eq — 0.160. The optimized field is shown 

in graph (a). The filter function ( ) /(w), scaled by 0.01, is plotted together with 

the spectrum in (b). The time-dependent occupation numbers (ground state ( ), 

first-excited state ( ) and third excited state ( — • — )) in (c) confirm that the 

transition occurs via the second excited state ( ). Due to the stronger constraints 

the convergence behaviour becomes oscillatory, shown in graph (d). 



is Eq = 0.400. The convergence of this optimization is shown in figure |10(b)[ After 
981 iterations, we obtain a target state occupation of 99.93%. The spectrum of the 
optimized pulse, shown in figure |10(a)[ exhibits contributions of all allowed frequency 
components. In particular, the zero- frequency component is dominant, being almost 
three times larger than the other frequency contributions. In the time-domain the zero- 
frequency component corresponds to a bias of e^v = —0.028 which is close to the value 
for which the potential becomes almost symmetric: e = —0.031. The optimized pulse 
together with e are shown in the middle panel of figure El 

The transfer process can be interpreted with the help of the following simplified 
picture: Assume the field would be almost static with e(t) ~ e. Then, the initial state 
becomes a superposition of the dressed states 

^{x,t = 0) = ^{vl-^l)e^'\ 

Small perturbations of the laser field around e rearrange the phases of the superposition 
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Figure 8. We apply the optimization algorithm with a double gaussian window, 
allowing only frequencies around wos and LO31 and in addition set Eq — 0.320. The 

optimized field is shown in graph (a). The filter function /(w) ( ), scaled by 0.01, 

is plotted together with the spectrum (b). The time-dependent occupation numbers 

(ground state ( ), first-excited state ( ) and second excited state ( )), 

shown in (c), confirm that the transition process is performed mainly via the third 
excited state ( — • — ). Due to the strong restrictions the convergence behaviour 
happens to be more oscillatory, shown in graph (d). 



SO that at the end of the pulse 

vl/(a;,T) = i=(^^ + ^De^'^. 

This superposition is located in the right well which completes the transfer. Note, that 
the phases ^1, Q2 are irrelevant in this case. 

The pulse we have obtained from the optimization is more difficult to explain since 
the oscillations around e are not small. To be able to analyse the transfer process 
in similar terms as discussed above, we have calculated the projections of the wave 
function \E'(x,t) onto the eigenfunctions of the Hamiltonian B.'" including the field e. 
These "dressed" occupation numbers are shown in the upper panel of figure El In 
the simplified interpretation we have implicitly assumed a complete localization in the 
left (right) well for the initial (target) state. Since this is not true for the potential 
chosen here, the dressed occupation numbers deviate slightly from 0.5, namely we have 
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|(^(x,0)|(^^)|2 = 0.57 and 0)|y?f)|2 = 0.43. At the end of the pulse we obtain 

the inverted occupation numbers, i.e. \{'^{x,T)\ipQ)\'^ = 0.45 and \{'^{x,T)\ipl)\'^ = 0.55 
which indicates the completed transfer (necessary condition). The laser pulse has also 
adjusted the phases of the expansion coefficients in the right way (sufficient condition): 
The relative phase difference of 0o 0i between t = and t = T was found to be 
0.8 * TT. This deviates from the simple picture where we would have expected a phase 
difference of vr. 

The transfer process described above has similarities with the one discovered in 
reference j2Hl where the authors have used an asymmetric double well to model a 
hydrogen transfer reaction. 




200 
t [a.u.] 



400 



Figure 9. In the lower panel we show the time-dependent occupation numbers of 

the ground-state ( ), the first excited state ( ) and the second excited state 

( ). The optimized laser field together with e — —0.031 is shown in the middle 

panel. In the top panel we plot the absolute values (squared) of the projections onto 
the two lowest eigenfunctions of the Hamiltonian including the field e. 



6. Conclusions 

We have presented a simple iterative scheme which allows for the optimization laser 
pulses under constraints on the spectrum of the laser and on its fluence. The scheme 
has been described in three different versions, one incorporating a given laser fluence, one 
restricting the spectrum of the laser pulse and the third one combining both constraints. 
Therefore, the scheme allows one to include realistic experimental constraints in the 
numerical optimization of laser pulses. 

To show that all three kinds of this scheme lead to high occupations of the target 
state we have applied them to drive the ground state of a ID asymmetric double well 
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Figure 10. Graph (a) shows the spectrum of the laser field ( ) and the rectangular 

filter ( ) scaled by 0.01. The convergence is shown in (b). The O indicates the 

iteration with the highest yield. 



potential to its first excited state. 

For all numerical tests we have obtained a high occupation > 99% in the target 
state. In the case of a fixed fluence {Eq = 0.080) we found a target state occupation of 
99.90% at the end of the pulse. Comparing the optimal laser pulses for different fiuences 
shows that, for fiuences larger than a certain critical value, target state occupations 
larger than 99% can always be achieved. With increasing fluence the optimized pulses 
employ a growing number of transition processes. Using spectral restrictions we are 
able to select between the different processes. We have calculated pulses that transfer 
the ground-state population to the target state only via the direct excitation process, 
explicitly without the direct process or via certain predefined intermediate levels. For 
the optimizations via an intermediate level we have additionally required a fixed laser 
fiuence. 

That it is possible to achieve a very high target state occupation with laser fields 
not containing any of the excitation frequencies has been clearly demonstrated by the 
last example. The laser spectrum was allowed to have frequency components only lower 
than the lowest resonance frequency. In addition we required the fiuence to be fixed. 
Like in the previous cases, the algorithm resulted in a laser pulse with a very high 
occupation of the target state. 

The results obtained in this work, clearly demonstrate that, in general, there exists 
no unique optimal laser pulse to achieve a given control target and that selection within 
the set of optimal pulses is possible by adding constraints to the optimization. With 
the methods presented here experimental constraints can be incorporated in the pulse 
optimization which makes the interpretation and analysis of the experimentally obtained 
laser pulses more reliable. The scheme allows one to study systematically the effects of 
different constraints on the target occupation and on the optimized laser field. Especially 
in the strong field regime, this leads to important insights in the various possible ways 
to achieve complete population transfer. 
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Appendix 

Appendix A. 1. Ineffectiveness of post- constraints 

To show the ineffectiveness of filtering the pulse spectrum after the optimization we have 
taken the optimized pulse from section IH. 1.11 and cut out the undesired frequencies. We 
then transform the pulse back to the time domain and propagate the time-dependent 
Schrodinger equation for the double well structure using this pulse. Since the pulse has 
a lower fluence after this procedure we also rescale the pulse so that 

/ dt t^{t) = 0.080. (A.l) 
Jo 

We apply this procedure to two cases: 

(i) Case 1: Restricition to the direct process. 

Filtering out all frequencies except u G [0.094, 0.236], which corresponds to consider 
only the direct process |0) — > results in a yield of 8% and after rescahng we 
obtain 44%. 

(ii) Case 2: Enforcing the indirect process. 

By filtering out all frequencies except u G [0.503, 0.833] we address only the indirect 
process |0) — * |2) Numerical propagation with the modified field results in 

a yield of 5% and 75% after rescaling. 

Comparing these numbers to the high yields found in sect ion 15.2.11 and 15.3. If we clearly 
see that filtering after optimization is ineffective. It is far more powerful to use the 
filtering in the optimization process. 



Appendix A. 2. Results from two-level system 

From the theory of two-level systems (or two-level atoms) [10] we can extract a good 
estimate for an optimal pulse, if the direct transition is allowed in dipole approximation. 
The estimate is extremely good, if no more than two-levels contribute to the process. 
This is the case if the excitation spectrum is well separated and the laser pulse is in 
the weak response regime. We have chosen a pulse length T = 400 which lies at the 
boundary of this regime but since the excitation energies are far apart from each other, 
we expect the two-level system to be a good approximation. The optimal pulse for 
a two- level-system (within the rotating wave-approximation(RWA)) that transfers all 
population from the ground-state to the excited state is a simple sinusoidal oscillation 

e{t) = Asm{uJoit) (A.2) 
where ujqi is the resonance frequency and A is the (optimal) amplitude given by: 

with the dipole matrix element /iqi = (0|/i|l) and T the length of the pulse. In our 
case we find A = 0.02003 and the corresponding fluence Eq = 0.0804. Applying this 
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pulse to the double well system, initially in the ground state, yields an occupation of 
99.30 % of the first excited state. In table IXTI we compare the results obtained from 
this simple estimate with the optimal control solution fixed to the same fluence. The 
results show that the two-level estimate is very successful for long times, however for 
short pulse lengths where T < 5 * 211/1001 it is not effective. This is due to the strength 
of the amplitude of the oscillation, it causes occupation also of the non-resonant levels. 

Table Al. Comparison of the yield P = |(l|^(T))p obtained with the two-level 
(RWA) pulse estimate versus the optimal control result. Note, that the period of the 
oscillation with the resonance frequency wqi is Tp — 40.08 a.u.. 
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